In [1]:
# import data
import geopandas as gpd
from shapely.geometry import Point, Polygon, MultiPoint
from matplotlib_scalebar.scalebar import ScaleBar
import matplotlib.pyplot as plt

import folium
from folium.plugins import MarkerCluster
import geopandas as gpd

import numpy as np

import gurobipy as gp
from gurobipy import GRB
In [2]:
# read data
stations = gpd.read_file("Data/Fire_Stations_SBC.shp")
county = gpd.read_file("Data/tl_2019_06083_faces.shp")
In [3]:
# remove water part
county =  county.loc[county["LWFLAG"] == 'L'].reset_index(drop=True)
county
Out[3]:
TFID STATEFP10 COUNTYFP10 TRACTCE10 BLKGRPCE10 BLOCKCE10 SUFFIX1CE ZCTA5CE10 UACE10 PUMACE10 ... METDIVFP CNECTAFP NECTAFP NCTADVFP LWFLAG OFFSET ATOTAL INTPTLAT INTPTLON geometry
0 257492722 06 083 002211 2 2004 None 93454 None 08301 ... None None None None L N 220098 +34.9497690 -120.3695124 POLYGON ((-120.37608 34.95466, -120.36430 34.9...
1 219296658 06 083 003102 1 1258 None 93436 None 08302 ... None None None None L N 255 +34.4938573 -120.4754525 POLYGON ((-120.47580 34.49391, -120.47569 34.4...
2 219296666 06 083 003102 1 1255 None 93436 None 08302 ... None None None None L N 1493 +34.4752880 -120.4521852 POLYGON ((-120.45261 34.47496, -120.45242 34.4...
3 219296667 06 083 003102 1 1257 None 93436 None 08302 ... None None None None L N 113683 +34.4706349 -120.4566533 POLYGON ((-120.45865 34.46769, -120.45850 34.4...
4 219296675 06 083 003102 1 1263 None 93436 None 08302 ... None None None None L N 166285 +34.4491085 -120.4587685 POLYGON ((-120.46231 34.45021, -120.46228 34.4...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
17274 219310979 06 083 001800 1 1089 None 93254 None 08302 ... None None None None L N 16549 +34.9985013 -119.8130570 POLYGON ((-119.81407 34.99936, -119.81388 34.9...
17275 219311368 06 083 001306 1 1014 None 93105 79282 08303 ... None None None None L N 149777 +34.4230712 -119.7456431 POLYGON ((-119.74828 34.42247, -119.74817 34.4...
17276 219313344 06 083 001906 4 4045 None 93105 None 08302 ... None None None None L N 175839 +34.5812649 -119.9848855 POLYGON ((-119.98708 34.58219, -119.98700 34.5...
17277 219313538 06 083 002906 2 2004 None 93117 79282 08303 ... None None None None L N 104502 +34.4500980 -119.8355449 POLYGON ((-119.83790 34.44967, -119.83785 34.4...
17278 219312610 06 083 002005 2 2001 None 93455 79417 08301 ... None None None None L N 115 +34.8613367 -120.4128687 POLYGON ((-120.41296 34.86140, -120.41282 34.8...

17279 rows × 45 columns

In [4]:
county_shape = county.unary_union
In [5]:
# plot data
A = stations.plot(color="red")
gpd.GeoSeries(county_shape).plot(ax=A)
Out[5]:
<Axes: >
No description has been provided for this image
In [6]:
stations.crs
Out[6]:
<Projected CRS: EPSG:2229>
Name: NAD83 / California zone 5 (ftUS)
Axis Info [cartesian]:
- X[east]: Easting (US survey foot)
- Y[north]: Northing (US survey foot)
Area of Use:
- name: United States (USA) - California - counties Kern; Los Angeles; San Bernardino; San Luis Obispo; Santa Barbara; Ventura.
- bounds: (-121.42, 32.76, -114.12, 35.81)
Coordinate Operation:
- name: SPCS83 California zone 5 (US Survey feet)
- method: Lambert Conic Conformal (2SP)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
In [7]:
county.crs
Out[7]:
<Geographic 2D CRS: EPSG:4269>
Name: NAD83
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: North America - onshore and offshore: Canada - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon. Puerto Rico. United States (USA) - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Hawaii; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands. British Virgin Islands.
- bounds: (167.65, 14.92, -40.73, 86.45)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich
In [8]:
# convert crs
county_shape = gpd.GeoDataFrame(geometry=[county.to_crs(stations.crs).unary_union], crs=stations.crs)
county_shape
Out[8]:
geometry
0 MULTIPOLYGON (((5776084.957 2033858.920, 57760...
In [9]:
# Now, let's plot them together.

A = stations.plot(color="red", zorder=100, markersize=10)
county_shape.plot(ax=A, facecolor='None', edgecolor="grey")

# Review for SCALEBAR!
scalebar = ScaleBar(dx=1, scale_formatter=lambda data, unit: f'{data},000 ft', length_fraction=0.25)
A.add_artist(scalebar)
Out[9]:
<matplotlib_scalebar.scalebar.ScaleBar at 0x21eac736410>
No description has been provided for this image
In [10]:
import folium
from folium.plugins import MousePosition
import geopandas as gpd


county_shape_4326 = county_shape.to_crs(4326)

# Calculate the centroid of the polygon
centroid = county_shape_4326.geometry.centroid.values[0].coords[0]

# Create a Folium map centered at the centroid of the polygon
m = folium.Map(location=[centroid[1], centroid[0]], zoom_start=9)

# Add the polygon to the map
folium.GeoJson(county_shape_4326, style_function=lambda x: { 'fillColor': 'transparent'}).add_to(m)

folium.GeoJson(stations.to_crs(4326)).add_to(m)

# Add MousePosition control to display cursor coordinates
#### Use the cursor coordinates to define your Area of Interest! ####
MousePosition().add_to(m)

# Display the map
m

C:\Users\18428\AppData\Local\Temp\ipykernel_13564\4168181903.py:9: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  centroid = county_shape_4326.geometry.centroid.values[0].coords[0]
Out[10]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [11]:
# changing boundary
max_x, min_x = -120.39153, -120.49805
max_y, min_y = 34.98936, 34.84762

# Create the rectangle geometry, which forms are study region boundary.

boundary = Polygon([(min_x, min_y), (min_x, max_y), (max_x, max_y), (max_x, min_y), (min_x, min_y)])
In [12]:
folium.GeoJson(boundary.__geo_interface__, style_function=lambda x: {'color': 'red', 'fillColor': 'transparent'}).add_to(m)

m
Out[12]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [13]:
boundary = gpd.GeoSeries(boundary, crs=4326).to_crs(stations.crs)[0]
boundary.area / 2.788e+7
Out[13]:
59.08324173317982
In [14]:
stations.within(boundary)
Out[14]:
0      False
1      False
2      False
3      False
4      False
       ...  
203    False
204    False
205    False
206    False
207    False
Length: 208, dtype: bool
In [15]:
study_stations =  stations.loc[stations.within(boundary)].reset_index(drop=True)
len(study_stations)
Out[15]:
9
In [16]:
boundary_county = county_shape.intersection(boundary)# the output will be assigned to boundary_county variable
boundary_county
Out[16]:
0    MULTIPOLYGON (((5817885.812 2191416.585, 58183...
dtype: geometry
In [36]:
COVERAGE_DISTANCE = 2 * 5280 # because 1*5280 feet is 1 mile 

# Draw buffer from each fire station in study_stations GeoDataFrame 
    # Buffer distance is COVERAGE_DISTANCE we've defined
stations_coverage = study_stations.buffer(COVERAGE_DISTANCE)

# Let's plot the county polgyon within our study region, and
    # stations_coverage buffers
A = stations_coverage.plot(facecolor="None")
boundary_county.plot(ax=A, zorder=-1)
plt.axis("off")

# How much area is covered?
    # To get to know, firstly take unary_union of the total station_coverage
total_stations_coverage = stations_coverage.unary_union
    # Then take the intersection with boundary_county, to exclude water area
total_stations_coverage = total_stations_coverage.intersection(boundary_county)
    # Then divide the area of the updated total_station_coverage with the boundary_county.area
coverage = total_stations_coverage.area / boundary_county.area
    # Multiply 100 to get the % value
coverage = coverage * 100

# Check out how much percentage of study region is covered by current fire stations.
    # Do they agree with each other? the coverage % and plot?
print(f'{round(coverage.values[0],2)} % covered')

76.7 % covered
No description has been provided for this image
In [37]:
redundant_coverage = 0
    # This two for-loops is the most efficient way to go through all 2 pair combinations of stations_coverage
for i in range(len(stations_coverage)-1):
    # Take one station coverage buffer
    cover_i = stations_coverage[i]
    
    for j in range(i+1, len(stations_coverage)):
        # Take another station coverage buffer which hasn't paired with i
        cover_j = stations_coverage[j]
        print("\t combination: ", i, j)
        
        # With the two stations combination, calculate the intersection area.
            # and sum up!
        redundant_coverage += cover_i.intersection(cover_j).area

print() # To add blank line
print("Redundant Coverage: ", redundant_coverage / 2.788e+7, " square miles")

	 combination:  0 1
	 combination:  0 2
	 combination:  0 3
	 combination:  0 4
	 combination:  0 5
	 combination:  0 6
	 combination:  0 7
	 combination:  0 8
	 combination:  1 2
	 combination:  1 3
	 combination:  1 4
	 combination:  1 5
	 combination:  1 6
	 combination:  1 7
	 combination:  1 8
	 combination:  2 3
	 combination:  2 4
	 combination:  2 5
	 combination:  2 6
	 combination:  2 7
	 combination:  2 8
	 combination:  3 4
	 combination:  3 5
	 combination:  3 6
	 combination:  3 7
	 combination:  3 8
	 combination:  4 5
	 combination:  4 6
	 combination:  4 7
	 combination:  4 8
	 combination:  5 6
	 combination:  5 7
	 combination:  5 8
	 combination:  6 7
	 combination:  6 8
	 combination:  7 8

Redundant Coverage:  75.52173527741837  square miles
In [38]:
def calculate_redundant_coverage(coverage_buffers):
    redundant_coverage = 0
        # This two for-loops is the most efficient way to go through all 2 pair combinations of stations_coverage
    for i in range(len(coverage_buffers)-1):
        # Take one station coverage buffer
        cover_i = coverage_buffers[i]

        for j in range(i+1, len(coverage_buffers)):
            # Take another station coverage buffer which hasn't paired with i
            cover_j = coverage_buffers[j]

            # With the two stations combination, calculate the intersection area.
                # and sum up!
            redundant_coverage += cover_i.intersection(cover_j).area

    print() # To add blank line
    print("Redundant Coverage: ", redundant_coverage / 2.788e+7, " square miles")
    return redundant_coverage / 2.788e+7
In [39]:
current_rc = calculate_redundant_coverage(stations_coverage)

Redundant Coverage:  75.52173527741837  square miles
In [40]:
boundary.bounds
Out[40]:
(5812247.644516614, 2139386.8749094545, 5845424.880169843, 2191733.8904521973)
In [41]:
minx, miny, maxx, maxy = boundary.bounds
print(minx, miny, maxx, maxy)

5812247.644516614 2139386.8749094545 5845424.880169843 2191733.8904521973
In [42]:
interval = 5280/3  # default: 5280/5 feet (1 mile / 5 = 0.2 mile) 


# Create arrays of x and y coordinates using np.arange
x_coords = np.arange(minx, maxx, interval)
y_coords = np.arange(miny, maxy, interval)

# Create a list to store the points
fishnet_points = []

# Generate points for the fishnet
for y in y_coords:
    for x in x_coords:
        fishnet_points.append((x, y))

# Print the number of points generated
print("Number of points in the fishnet:", len(fishnet_points))

fishnet_points = gpd.GeoSeries([Point(pt_cd) for pt_cd in fishnet_points])
fishnet_points.plot(figsize=(15,5))

Number of points in the fishnet: 570
Out[42]:
<Axes: >
No description has been provided for this image
In [43]:
fishnet_points = fishnet_points.loc[fishnet_points.intersects(boundary_county[0])]
fishnet_points = fishnet_points.loc[fishnet_points.intersects(total_stations_coverage[0])].reset_index(drop=True)

fishnet_points.plot(figsize=(15,5))
Out[43]:
<Axes: >
No description has been provided for this image
In [44]:
A = fishnet_points.plot(figsize=(15,5), color="black", markersize=10)
study_stations.plot(ax=A, color="blue", marker="*", markersize=400)
Out[44]:
<Axes: >
No description has been provided for this image
In [45]:
# optimize the problem
num_demands = len(fishnet_points) 
num_facilities = len(study_stations)

D_ij = [[np.sqrt((fishnet_points.x[i] - study_stations.geometry.x[j])**2 + (fishnet_points.y[i] - study_stations.geometry.y[j])**2) 
         for j in range(num_facilities)] for i in range(num_demands)]

N_i = dict(zip(range(num_demands), [[] for _ in range(num_demands)]))

for i in range(num_demands):
    for j in range(num_facilities):
        if D_ij[i][j] <= COVERAGE_DISTANCE:
            N_i[i].append(j)

# Create a new model
model = gp.Model("LSCP")

# Decision variables
# x[i] = 1 if facility i is selected, 0 otherwise
x = model.addVars(num_facilities, vtype=GRB.BINARY, name="x")

# Objective function: minimize the number of selected facilities
model.setObjective(gp.quicksum(x[j] for j in range(num_facilities)), GRB.MINIMIZE)

# Constraints
# Each demand point must be covered by at least one selected facility
for i in range(num_demands):
    model.addConstr(gp.quicksum(x[j] for j in N_i[i]) >= 1, name=f"cover_demand_point_{j}")

# Optimize the model
model.optimize()

Gurobi Optimizer version 11.0.1 build v11.0.1rc0 (win64 - Windows 11.0 (22631.2))

CPU model: 11th Gen Intel(R) Core(TM) i7-11370H @ 3.30GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 383 rows, 9 columns and 794 nonzeros
Model fingerprint: 0xbcb456a9
Variable types: 0 continuous, 9 integer (9 binary)
Coefficient statistics:
  Matrix range     [1e+00, 1e+00]
  Objective range  [1e+00, 1e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]
Found heuristic solution: objective 8.0000000
Presolve removed 383 rows and 9 columns
Presolve time: 0.00s
Presolve: All rows and columns removed

Explored 0 nodes (0 simplex iterations) in 0.01 seconds (0.00 work units)
Thread count was 1 (of 8 available processors)

Solution count 1: 8 

Optimal solution found (tolerance 1.00e-04)
Best objective 8.000000000000e+00, best bound 8.000000000000e+00, gap 0.0000%
In [46]:
selected_bool = [x[i].x > 0.5 for i in range(num_facilities)]
selected_facilities = study_stations.loc[selected_bool].reset_index(drop=True)
unselected_facilities = study_stations.loc[~np.array(selected_bool)].reset_index(drop=True)

print("Objective value: ", model.objVal, sum(selected_bool))
print("Current running stations: ", len(study_stations))

Objective value:  8.0 8
Current running stations:  9
In [47]:
A = gpd.GeoSeries(MultiPoint(fishnet_points)).plot(figsize=(15,5), color="grey", markersize=5)
unselected_facilities.plot(ax=A, color="blue", legend=True, label="Unselected", marker="*", markersize=300, zorder=10)
selected_facilities.plot(ax=A, color="red", marker="*", markersize=300, legend=True, label="Selected")
selected_facilities.buffer(COVERAGE_DISTANCE).plot(ax=A, edgecolor="red", facecolor="None", linestyle="dashed")
plt.legend()
plt.axis("off")
Out[47]:
(5815990.983205412, 5853511.824250226, 2132101.015651713, 2202523.549244632)
No description has been provided for this image
In [48]:
# Convert GeoPandas GeoSeries to GeoDataFrame
fishnet_gdf = gpd.GeoDataFrame(geometry=fishnet_points, crs=stations.crs).to_crs(4326)

# Create a base map
m = folium.Map(location=[fishnet_gdf.unary_union.centroid.y, fishnet_gdf.unary_union.centroid.x], zoom_start=11)



# Plot fishnet
for idx, row in fishnet_gdf.iterrows():
    folium.CircleMarker(location=[row.geometry.y, row.geometry.x],
                        radius=1,
                        color='grey', fillcolor="grey").add_to(m)

# Plot study stations
for idx, row in unselected_facilities.to_crs(4326).iterrows():
    folium.Marker(location=[row.geometry.y, row.geometry.x],
                   icon=folium.Icon(color='blue'),
                   popup=row['Label']).add_to(m)

# Plot selected facilities
for idx, row in selected_facilities.to_crs(4326).iterrows():
    folium.Marker(location=[row.geometry.y, row.geometry.x],
                   icon=folium.Icon(color='red', icon='star'),
                   popup=row['Label']).add_to(m)

# Plot buffer around selected facilities
for idx, row in selected_facilities.to_crs(4326).iterrows():
    folium.vector_layers.Circle(location=[row.geometry.y, row.geometry.x],
                                 radius=COVERAGE_DISTANCE/3.281, # feet to meter conversion
                                 color='red',
                                 fill=False,
                                 dash_array='10').add_to(m)

# Display the map
m
Out[48]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [49]:
optimized_rc =  calculate_redundant_coverage(selected_facilities.buffer(COVERAGE_DISTANCE))

Redundant Coverage:  51.62504510669677  square miles
In [50]:
spatial_efficiency = model.objVal / len(study_stations) * 100

print("Spatial Efficienty is ", round(spatial_efficiency, 2), " %")

Spatial Efficienty is  88.89  %
In [51]:
print("Coverage Distance (mi): " , COVERAGE_DISTANCE) # Double Check this value before putting into table ** 
print("Current Stations #: " ,len(study_stations))
print("Optimized Stations #: ", len(selected_facilities))
print("Current Redundant Coverage (sqmi): ", current_rc)
print("Optimized Redundant Coverage (sqmi): ", optimized_rc)
print("Spatial Efficiency (%): ", spatial_efficiency)

Coverage Distance (mi):  10560
Current Stations #:  9
Optimized Stations #:  8
Current Redundant Coverage (sqmi):  75.52173527741837
Optimized Redundant Coverage (sqmi):  51.62504510669677
Spatial Efficiency (%):  88.88888888888889
In [ ]: